function [B_ols, T_stat, BIC, RSS, Fore] = var_est_f(x , lags, c ,tl, h ,opt)

% OLS estimates
[T N] = size(x);
Y = reshape(x(lags+1:end,:),(T-lags)*N,1);
xx = [];
for j = 1:lags;
    xx = [xx x(lags+1-j:end-j,:) ];
end;

if c == 1;
    xx = [ones(T-lags,1) xx];
end
X = kron(eye(N), xx);

B_ols = inv(X'*X)*X'*Y;
res = Y - X*B_ols;
RES=reshape(res,(T-lags)*N/N,N);
RSS = RES'*RES;
c_cov = RES'*RES/(size(RES,1)-lags-c);
se = sqrt(diag(inv(X'*X)*kron(diag(diag(c_cov)),eye(size(B_ols,1)/N))));
T_stat = B_ols./se;

BIC = log(det(c_cov))+(N+N^2*(lags))*log(length(RES))/length(RES);

% Companion matrix
AA = reshape(B_ols,size(B_ols,1)/N,N)';
A = [AA(:,c+1:end); [eye(N*(lags-1)) zeros(N*(lags-1),N)]];


if c~0;
    CC = zeros(lags*N,1);
    for jj = 1:h
        CC = CC + A^(jj-1)*[B_ols(1:lags*N+1:end,1) ;zeros(N*(lags-1),1)];
    end

else
    CC = zeros(lags*N,1);
end


fo = CC + A^h*[x(end,1:N) xx(end,c+1:end-N)]';

Fore = fo(1:N,1);


if opt ==1
    % IRF
    M = zeros((lags)); M(1,1)=1;
    CC = kron(M,chol(c_cov)');
    for i=0:tl-1
        irf(:,:,i+1)=A^i*CC;
    end

    % Figure

    labels{1} = 'INFL';
    labels{2} = 'GDP';
    labels{3} = 'FFR';

    figure(1);
    count = 0 ;
    for i = 1:N
        for j = 1:N
            count = count +1;
            subplot(N,N,count); plot((1:tl)',squeeze(irf(i,j,:)),'LineWidth',2);grid on ; axis([1 tl -inf inf])
            title([labels{i},' to ',labels{j}],'Color',[0 0 1],'FontSize',11)
        end
    end

set(gcf,'Color',[0.9725 0.9725 0.9725])


count = 0;
a=cumsum(irf.^2,3);
figure(2);
for i = 1:N
    v = squeeze(a(:,i,:))./squeeze(sum(a,2));
    for j=1:N
        count = count +1;
        subplot(N,N,i+N*(j-1)); plot(v(j,:)*100,'LineWidth',2);axis([-inf inf 0 100 ])
        title([labels{i},' to ',labels{j}],'Color',[0 0 1],'FontSize',11)
    end

end
set(gcf,'Color',[0.9725 0.9725 0.9725])
end